home *** CD-ROM | disk | FTP | other *** search
/ Technotools / Technotools (Chestnut CD-ROM)(1993).ISO / lang_oth / linpklib / schex.for < prev    next >
Text File  |  1985-01-12  |  7KB  |  248 lines
  1.       SUBROUTINE SCHEX(R,LDR,P,K,L,Z,LDZ,NZ,C,S,JOB)
  2.       INTEGER LDR,P,K,L,LDZ,NZ,JOB
  3.       REAL R(LDR,1),Z(LDZ,1),S(1)
  4.       REAL C(1)
  5. C
  6. C     SCHEX UPDATES THE CHOLESKY FACTORIZATION
  7. C
  8. C                   A = TRANS(R)*R
  9. C
  10. C     OF A POSITIVE DEFINITE MATRIX A OF ORDER P UNDER DIAGONAL
  11. C     PERMUTATIONS OF THE FORM
  12. C
  13. C                   TRANS(E)*A*E
  14. C
  15. C     WHERE E IS A PERMUTATION MATRIX.  SPECIFICALLY, GIVEN
  16. C     AN UPPER TRIANGULAR MATRIX R AND A PERMUTATION MATRIX
  17. C     E (WHICH IS SPECIFIED BY K, L, AND JOB), SCHEX DETERMINES
  18. C     A ORTHOGONAL MATRIX U SUCH THAT
  19. C
  20. C                           U*R*E = RR,
  21. C
  22. C     WHERE RR IS UPPER TRIANGULAR.  AT THE USERS OPTION, THE
  23. C     TRANSFORMATION U WILL BE MULTIPLIED INTO THE ARRAY Z.
  24. C     IF A = TRANS(X)*X, SO THAT R IS THE TRIANGULAR PART OF THE
  25. C     QR FACTORIZATION OF X, THEN RR IS THE TRIANGULAR PART OF THE
  26. C     QR FACTORIZATION OF X*E, I.E. X WITH ITS COLUMNS PERMUTED.
  27. C     FOR A LESS TERSE DESCRIPTION OF WHAT SCHEX DOES AND HOW
  28. C     IT MAY BE APPLIED, SEE THE LINPACK GUIDE.
  29. C
  30. C     THE MATRIX Q IS DETERMINED AS THE PRODUCT U(L-K)*...*U(1)
  31. C     OF PLANE ROTATIONS OF THE FORM
  32. C
  33. C                           (    C(I)       S(I) )
  34. C                           (                    ) ,
  35. C                           (    -S(I)      C(I) )
  36. C
  37. C     WHERE C(I) IS REAL, THE ROWS THESE ROTATIONS OPERATE ON
  38. C     ARE DESCRIBED BELOW.
  39. C
  40. C     THERE ARE TWO TYPES OF PERMUTATIONS, WHICH ARE DETERMINED
  41. C     BY THE VALUE OF JOB.
  42. C
  43. C     1. RIGHT CIRCULAR SHIFT (JOB = 1).
  44. C
  45. C         THE COLUMNS ARE REARRANGED IN THE FOLLOWING ORDER.
  46. C
  47. C                1,...,K-1,L,K,K+1,...,L-1,L+1,...,P.
  48. C
  49. C         U IS THE PRODUCT OF L-K ROTATIONS U(I), WHERE U(I)
  50. C         ACTS IN THE (L-I,L-I+1)-PLANE.
  51. C
  52. C     2. LEFT CIRCULAR SHIFT (JOB = 2).
  53. C         THE COLUMNS ARE REARRANGED IN THE FOLLOWING ORDER
  54. C
  55. C                1,...,K-1,K+1,K+2,...,L,K,L+1,...,P.
  56. C
  57. C         U IS THE PRODUCT OF L-K ROTATIONS U(I), WHERE U(I)
  58. C         ACTS IN THE (K+I-1,K+I)-PLANE.
  59. C
  60. C     ON ENTRY
  61. C
  62. C         R      REAL(LDR,P), WHERE LDR.GE.P.
  63. C                R CONTAINS THE UPPER TRIANGULAR FACTOR
  64. C                THAT IS TO BE UPDATED.  ELEMENTS OF R
  65. C                BELOW THE DIAGONAL ARE NOT REFERENCED.
  66. C
  67. C         LDR    INTEGER.
  68. C                LDR IS THE LEADING DIMENSION OF THE ARRAY R.
  69. C
  70. C         P      INTEGER.
  71. C                P IS THE ORDER OF THE MATRIX R.
  72. C
  73. C         K      INTEGER.
  74. C                K IS THE FIRST COLUMN TO BE PERMUTED.
  75. C
  76. C         L      INTEGER.
  77. C                L IS THE LAST COLUMN TO BE PERMUTED.
  78. C                L MUST BE STRICTLY GREATER THAN K.
  79. C
  80. C         Z      REAL(LDZ,NZ), WHERE LDZ.GE.P.
  81. C                Z IS AN ARRAY OF NZ P-VECTORS INTO WHICH THE
  82. C                TRANSFORMATION U IS MULTIPLIED.  Z IS
  83. C                NOT REFERENCED IF NZ = 0.
  84. C
  85. C         LDZ    INTEGER.
  86. C                LDZ IS THE LEADING DIMENSION OF THE ARRAY Z.
  87. C
  88. C         NZ     INTEGER.
  89. C                NZ IS THE NUMBER OF COLUMNS OF THE MATRIX Z.
  90. C
  91. C         JOB    INTEGER.
  92. C                JOB DETERMINES THE TYPE OF PERMUTATION.
  93. C                       JOB = 1  RIGHT CIRCULAR SHIFT.
  94. C                       JOB = 2  LEFT CIRCULAR SHIFT.
  95. C
  96. C     ON RETURN
  97. C
  98. C         R      CONTAINS THE UPDATED FACTOR.
  99. C
  100. C         Z      CONTAINS THE UPDATED MATRIX Z.
  101. C
  102. C         C      REAL(P).
  103. C                C CONTAINS THE COSINES OF THE TRANSFORMING ROTATIONS.
  104. C
  105. C         S      REAL(P).
  106. C                S CONTAINS THE SINES OF THE TRANSFORMING ROTATIONS.
  107. C
  108. C     LINPACK. THIS VERSION DATED 08/14/78 .
  109. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  110. C
  111. C     SCHEX USES THE FOLLOWING FUNCTIONS AND SUBROUTINES.
  112. C
  113. C     BLAS SROTG
  114. C     FORTRAN MIN0
  115. C
  116.       INTEGER I,II,IL,IU,J,JJ,KM1,KP1,LMK,LM1
  117.       REAL RJP1J,T
  118. C
  119. C     INITIALIZE
  120. C
  121.       KM1 = K - 1
  122.       KP1 = K + 1
  123.       LMK = L - K
  124.       LM1 = L - 1
  125. C
  126. C     PERFORM THE APPROPRIATE TASK.
  127. C
  128.       GO TO (10,130), JOB
  129. C
  130. C     RIGHT CIRCULAR SHIFT.
  131. C
  132.    10 CONTINUE
  133. C
  134. C        REORDER THE COLUMNS.
  135. C
  136.          DO 20 I = 1, L
  137.             II = L - I + 1
  138.             S(I) = R(II,L)
  139.    20    CONTINUE
  140.          DO 40 JJ = K, LM1
  141.             J = LM1 - JJ + K
  142.             DO 30 I = 1, J
  143.                R(I,J+1) = R(I,J)
  144.    30       CONTINUE
  145.             R(J+1,J+1) = 0.0E0
  146.    40    CONTINUE
  147.          IF (K .EQ. 1) GO TO 60
  148.             DO 50 I = 1, KM1
  149.                II = L - I + 1
  150.                R(I,K) = S(II)
  151.    50       CONTINUE
  152.    60    CONTINUE
  153. C
  154. C        CALCULATE THE ROTATIONS.
  155. C
  156.          T = S(1)
  157.          DO 70 I = 1, LMK
  158.             CALL SROTG(S(I+1),T,C(I),S(I))
  159.             T = S(I+1)
  160.    70    CONTINUE
  161.          R(K,K) = T
  162.          DO 90 J = KP1, P
  163.             IL = MAX0(1,L-J+1)
  164.             DO 80 II = IL, LMK
  165.                I = L - II
  166.                T = C(II)*R(I,J) + S(II)*R(I+1,J)
  167.                R(I+1,J) = C(II)*R(I+1,J) - S(II)*R(I,J)
  168.                R(I,J) = T
  169.    80       CONTINUE
  170.    90    CONTINUE
  171. C
  172. C        IF REQUIRED, APPLY THE TRANSFORMATIONS TO Z.
  173. C
  174.          IF (NZ .LT. 1) GO TO 120
  175.          DO 110 J = 1, NZ
  176.             DO 100 II = 1, LMK
  177.                I = L - II
  178.                T = C(II)*Z(I,J) + S(II)*Z(I+1,J)
  179.                Z(I+1,J) = C(II)*Z(I+1,J) - S(II)*Z(I,J)
  180.                Z(I,J) = T
  181.   100       CONTINUE
  182.   110    CONTINUE
  183.   120    CONTINUE
  184.       GO TO 260
  185. C
  186. C     LEFT CIRCULAR SHIFT
  187. C
  188.   130 CONTINUE
  189. C
  190. C        REORDER THE COLUMNS
  191. C
  192.          DO 140 I = 1, K
  193.             II = LMK + I
  194.             S(II) = R(I,K)
  195.   140    CONTINUE
  196.          DO 160 J = K, LM1
  197.             DO 150 I = 1, J
  198.                R(I,J) = R(I,J+1)
  199.   150       CONTINUE
  200.             JJ = J - KM1
  201.             S(JJ) = R(J+1,J+1)
  202.   160    CONTINUE
  203.          DO 170 I = 1, K
  204.             II = LMK + I
  205.             R(I,L) = S(II)
  206.   170    CONTINUE
  207.          DO 180 I = KP1, L
  208.             R(I,L) = 0.0E0
  209.   180    CONTINUE
  210. C
  211. C        REDUCTION LOOP.
  212. C
  213.          DO 220 J = K, P
  214.             IF (J .EQ. K) GO TO 200
  215. C
  216. C              APPLY THE ROTATIONS.
  217. C
  218.                IU = MIN0(J-1,L-1)
  219.                DO 190 I = K, IU
  220.                   II = I - K + 1
  221.                   T = C(II)*R(I,J) + S(II)*R(I+1,J)
  222.                   R(I+1,J) = C(II)*R(I+1,J) - S(II)*R(I,J)
  223.                   R(I,J) = T
  224.   190          CONTINUE
  225.   200       CONTINUE
  226.             IF (J .GE. L) GO TO 210
  227.                JJ = J - K + 1
  228.                T = S(JJ)
  229.                CALL SROTG(R(J,J),T,C(JJ),S(JJ))
  230.   210       CONTINUE
  231.   220    CONTINUE
  232. C
  233. C        APPLY THE ROTATIONS TO Z.
  234. C
  235.          IF (NZ .LT. 1) GO TO 250
  236.          DO 240 J = 1, NZ
  237.             DO 230 I = K, LM1
  238.                II = I - KM1
  239.                T = C(II)*Z(I,J) + S(II)*Z(I+1,J)
  240.                Z(I+1,J) = C(II)*Z(I+1,J) - S(II)*Z(I,J)
  241.                Z(I,J) = T
  242.   230       CONTINUE
  243.   240    CONTINUE
  244.   250    CONTINUE
  245.   260 CONTINUE
  246.       RETURN
  247.       END
  248.